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Abstract 


Large station displacements observed from Imperial Valley GPS campaigns are attributed to the 
November 24, 1987 Superstition Hills earthquake sequence. Thirty sites from a 42 station GPS network 
established in 1986 have been reoccupied during 1988 and/or 1990. Displacements at three sites within 3 
kilometers of the surface rupture approach 0.5 m. Eight additional stations within 20 km of the seismic 
zone are displaced at least 10 cm. This is the first occurrence of a large earthquake (Mg 6.6) within a 
pr eexi sting GPS network. Best-fitting uniform slip models rectangular dislocations in an elastic half- 
space indicate 130 ± 8 cm right-lateral displacement along the northwest-trending Superstition Hills fault 
and 30 ± 10 cm left-lateral displacement along the conjugate northeast-trending Elmore Ranch fault. The 
geodetic moments are 9.4 x 10 25 dyne-cm and 2.3 x 10 25 dyne- cm for the Superstition Hills and Elmore 
Ranch faults, respectively, consistent with teleseismic source parameters. The data also suggest that 
postseismic slip along the Superstition Hills fault is concentrated at shallow depths. Distributed slip 
solutions using Singular Value Decomposition indicate near uniform displacement along the Elmore Ranch 
fault and concentrated slip to the northwest and southeast along the Superstition Hills fault. A significant 
component of non-seismic secular displacement is observed across the Imperial Valley, which is attributed 
in part to interseismic plate- boundary deformation. 
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1. Introduction 

The Global Positioning System (GPS) is rapidly becoming one of the most important 
tools for studying tectonic deformation [e.g., Larsen, 1001]. By recording signals from 
earth orbiting satellites It is possible to determine 3-dimensional coordinates of geodetic 
monuments ■with high accuracy [e.g., King et al., 1085; Wells et al., 1087]. The station 
displacement or deformation between surveys is measured when observations are 
repeated. GPS can be used to monitor secular deformation such as that associated with 
plate motion, or to record rapid strain fluctuations such as those due to seismic and 
volcanic activity. In its final configuration scheduled for the mid 1990’s, 21 satellites will 
be deployed in 6 orbital planes (with 3 additional satellites used as active spares). When 
G3PS becomes fully operational, it will be possible to continuously determine 3- 
dimensional positions anywhere on or near the earth. The available satellite 
constellation existing for the last several years was optimized for North America, making 
geodetic studies in California practical. The observation window in which enough 
satellites have been visible to obtain the high accuracies necessary to measure tectonic 
motion has been about 6 to 8 hours each day. 

On November 24, 1987, two large earthquakes separated by 12 hours occurred in the 
northwest portion of the Imperial Valley region of southern California. The first event 
was located on a northeast-trending seismic lineament and was followed 12 hours by 
rupture along the northwest-trending Superstition hills fault. What makes this 
earthquake sequence so significant from a GPS standpoint, is that it occurred spatially 
and temporally within a preexisting GPS network. This network was established in the 
Imperial Valley in 1986, with partial resurveys in 1988 and 1990. Fifteen stations are 
located within 20 km of the rupture zone; three stations are within 3 km. This is the 
first occurrence of a large earthquake within a preexisting GPS network. 
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We compute GPS displacements in the Imperial Valley between 1986 and 1990. 
Observed station movements of nearly 0.5 meters are attributed to the Superstition Hills 
earthquake sequence. The earthquake-induced displacements are inverted to estimate 
seismic slip and the geodetic moment along the rupture planes. In addition, there is a 
large component of deformation which can not be explained by the seismic disturbance. 
We assume this to be, in part, a manifestation of continuous strain accumulation across 
the Imperial Valley due to the relative motion of the Pacific and North American plates. 

2. Imperial Valley Seismicity and Tectonics 

The Imperial Valley region of southern California is a complex transition zone 
between crustal spreading in the Gulf of California and rightrlateral transform motion 
along the San Andreas fault (Figure 1) [Lomnitz et al, 1970; Elders et al. , 1972]. The 
valley is 4 - 12 million years old and has been filled by up to 15 km of late Cenozoic 
sediments [ Larson et al., 1968; Moore and Buffington, 1968; Ingle, 1974; Fuis et al., 1982]. 
The major fault systems and structural grain of the valley trend to the northwest, 
roughly parallel to the direction of plate motion. A significant fraction of the North 
American and Pacific relative motion may be accommodated across this region. 

The Imperial Valley is one of the most seismically active regions of California (Figure 
1). Much of this activity occurs along the Imperial fault and within the Brawley Seismic 
Zone [Johnson and Hill, 1982]. Several large earthquakes have occurred in and near the 
Imperial Valley since 1940 (Figure 1). The Imperial fault ruptured with a Mg 7.1 event 
in 1940 and a Ml 6.6 event in 1979 (17.5. Geol. Surv., 1982]. Segments of the San 
Jacinto fault system broke with a Ml 6.2 earthquake in 1954 and a Ml 6.5 event in 1968 
(Borrego Mountain). The most GPS relevant episode of seismic activity occurred in 1987 
along the Superstition Hills segment of the San Jacinto fault system, with a Mg 6.2 


ear thquake on a northeast-trending seismic lineament followed 12 hours by a M s 6.6 
event on the Superstition Hills fault. 

Conventional geodetic measurements suggest considerable deformation across the 
Imperial Valley. In fact, a significant fraction of the Pacific-North American relative 
plate motion may be accommodated here. New global plate model estimates (NUVEL-1) 
[DeMets et aL , 1987; DeMets et al., 1990] predict the rate of motion between the North 
American and Pacific plates (at Imperial Valley coordinates: 33.0 *N, 115.5 *W) is 4.7 
cm/yr oriented N39.6*W. Triangulation measurements spanning this region suggest 4.3 
cm/yr of plate-boundary deformation averaged between 1941 and 1986 [Snay and Drew , 
1988]. Trilateration measurements made by the U.S.G.S. between 1972 and 1987 
indicate 3.45 cm/yr relative movement between stations on opposite sides of the valley 
[Prescott et al., 1987b]. The orientations of the conventional geodetic displacements are 
approximately N40 * W, although to some extent the direction is non-unique and depends 
upon a priori assumptions [e.g., Prescott, 1981]. In addition, the conventional geodetic 
measurements indicate that deformation is concentrated in a narrow 20 km wide zone 
along the Imperial fault, while to the north it is distributed over a region at least 50 km 
wide. Presumably, deformation is transferred from the Imperial fault, which acts as the 
primary strain release mechanism near the border, to distributed shear along the San 
Andreas, San Jacinto, and Elsinore faults. 

3. Superstition Hills Earthquake Sequence 

. m 

On November 24, 1987 (1:54 GMT), a M s 6.2 earthquake occurred along a 
northeast-trending seismic lineament north east of t he Superstition Hills fault (Figure 2) 
[Magistrate et al, 1988]. The focal mechanism and aftershock sequence, which extended 
for 26 km into the Brawley Seismic Zone, are consistent with left-lateral strike slip 



motion on a vertical fault. Seven foreshocks were recorded in the 22 minutes prior to the 
main event, including two with M L > 4.0. Surface rupture consisted of a complex 
pattern of left-lateral northeast-trending offsets ranging in length from 1.5 to 10 km, and 
with maximum displacements between 3 and 13 cm [Budding and Sharp, 1988; Hudnut et 
al. , 1989a]. We refer to this northeast-trending lineament as the Elmore Ranch fault, 
although more precisely this name refers only to the longest of the surface fractures. 

Twelve hours after the Elmore Ranch event (13:15 GMT), a M s 6.6 earthquake 
occurred along the northwest-trending Superstition Hills fault. The epicenter was near 
the intersection of the Elmore Ranch and Superstition Hills faults. Strong ground 
motion and teleseismic data suggest the rupture process for this second event was 
somewhat complicated, consisting of multiple sub-events [Bent et al, 1989; Frankel and 
Wennerberg, 1989; Hwang et al., 1990; Wald et al., 1990]. Surface rupture extended 24 
lfTw along the previously mapped trace of the fault [ Williams and Magistrate, 1989]. Up 
to 50 cm right-lateral displacement were measured immediately following the earthquake. 
The aftershock pattern was concentrated slightly to the west of the fault and did not 
extend the length of the surface rupture. Magistrate et al. [1989] suggested the aftershock 
sequence was highly correlated with basement structure. Both the Superstition Hills and 
Elmore Ranch events triggered sympathetic surface offsets along the Imperial, San 
Andreas, and Coyote Creek faults [McGill et al., 1989; Hudnut and Clark, 1989]. 

Significant afterslip was recorded along the Superstition Hills fault following the 
second mainshock [ Williams and Magistrate, 1989]. No afterslip was measured along any 
of the surface ruptures associated with the Elmore Ranch event. In fact, all seismic 
activity essentially stopped along this segment after the initiation of the 2nd mainshock. 

One of the most interesting aspects of this earthquake sequence is the conjugate 
nature of faulting. That is, two surface ruptures oriented nearly perpendicular to each 
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other. As discussed below, this type of fault interaction may be typical of Imperial 
Valley tectonics and may dictate the mode of stress/strain transfer from one fault system 
to another. 

Wha t makes the Superstition Hills earthquake sequence unique from a GPS 
perspective is that it occurred within a preexisting GPS network. Stations located near 
the seismic rupture zone were displaced nearly 0.5 meters. These movements, as well as 
smaller displacements observed at nearby stations, can be Inverted to infer properties of 
the rupture process. This is the first time GPS measurements have recorded the 
deformation from a large earthquake. 

4. GPS Observations 

The signal structure broadcast from each GPS satellite consists of two carrier phase 
signals modulated by a navigational message and pseudo-random codes. The two carrier 
frequencies, known as the LI and L2 phases, are broadcast at 1575.42 Mhz (Ll) and 
1227.60 Mhz (L2). This is equivalent to wavelengths of about 19 cm for the Ll and 24 
cm for the L2. The navigational message contains the satellite coordinates (broadcast 
ephemeris), clock parameters, satellite health, and general system status. The pseudo- 
random codes are accurate time marks which allow a GPS receiver to determine the 
transmission time of the signal. When scaled by the speed of light, the pseudorange, or 
the time-biased satellite-receiver distance is computed. If measurements to at least 4 
GPS satellites are available, and if the satellite coordinates are known (usually with the 
broadcast ephemeris), the 3-dimensional receiver position and the satellite-receiver time 
offset can be determined. The positioning accuracy with the pseudorange is 1 to 100 m, 
depending on whether the P or C/A code Is used, the receiver type, length of 
observation, and the static or kinematic behavior of the instrument. It is the 


pseudorange which will be used for civilian and military navigation. For highly accurate 
geodetic positioning, however, it is necessary to use the carrier phase measurements in a 
post-processing mode. That is, the data collected in the field are brought back to the 
office or laboratory for analysis, usually with a fairly robust computer software system. 
More complete details about the Global Positioning System, including theoretical aspects 
and processing methods, are found in King et al. [1085], Wells et at. [1987], and Rocken 
[1088]. 

GPS Surveys — Data Collection 

The GPS data for this study were collected during 4 Imperial Valley field campaigns 
from 1986 to 1900 (Table 1). A total of 46 stations in or near the valley have been 
occupied at least once during this interval, while 30 sites have been reoccupied since 1986 
(Figure 3). TI-4100 GPS receivers supporting GESAR software were used during 1986 
and 1988, while Trimble 4000SST instruments were used during the 1990 survey. A 
comprehensive discussion regarding the data collection methods is given in Larsen [1991]. 

The National Geodetic Survey (NGS) occupied 54 sites in southern California in 
May/June 1986; 42 stations were located in or near the Imperial Valley. Most marks 
were observed at least two days, although redundant baselines were relatively uncommon 
(i.e., simultaneous occupation of the same station-station pair for two or more days). 
This was one of the first GPS networks established to investigate crustal motions. In 
fact, GPS surveying at this time was still at an "experimental stage." The scheduled 4.5 
hour daily occupations were somewhat less than the 6-8 hour sessions typical of more 
recent southern California campaigns. There was less than a 2 hour period each day in 
which more than 3 satellites were simultaneously tracked, and there was a scheduled 1 
hour gap in Satellite 6 during the middle of the measurement session. In addition, due 
to a variety of equipment and logistical problems, less than 2 hours of data were 
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collected at several sites. The data were generally noisy and processing was unusually 
tedious. The final double difference phase observables may contain uncorrected cycle 
slips. Because of these unfortunate circumstances due to the earliness of this survey, the 
quality of the data is fairly poor. 

In February /March 1988 university field crews (UNAVCO) observing for 9 days 
reoccupied 15 of the Imperial Valley marks, as well as establishing 4 new stations along 
the Salton Sea. These new sites were installed at tide gauge monitoring facilities and 
will be used to constrain vertical GPS baseline components and for monitoring tectonic 
tilting in the Salton Trough. Most monuments were occupied for 2 days (1 was observed 
3 days, 1 for 4 days, and 2 for 1 day). The scheduled nightly observation scenario lasted 
7.5 hours, with a total of 7 satellites tracked. 

In March/April 1988 the NGS reoccupied 21 of the 1986 stations (7 of which were 
observed by the university crews a month earlier). Most sites were occupied only 1 day. 
The daily observation period of 6.0 hours was slightly shorter than the February /March 
survey, although 7 satellites were tracked each day. 

In addition to the UNAVCO and NGS campaigns, GPS observations in 1988 were 
made at Mojave (California), Westford (Massachusetts), and Richmond (Florida). These 
sites are continuously monitored as part of the Cooperative International GPS Network 
(CIGNET) {Chin, 1988]. Data from CIGNET stations are used to improve satellite 
orbits in the GPS processing. Unfortunately, these observations are frequently of poor 
quality and contain Abundant cycle slips. During the 1988 campaign, observations were 
not always available or usable at all CIGNET sites on all days {Larsen, 1991]. 

Station MOUN (Mound), which was surveyed in 1986, was destroyed from the 1987 
earthquake sequence. The site is located less than 1 km from the surface rupture of the 
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Superetition Hills fault. Field investigation during early 1088 reveiled that the 
monument and supporting concrete base had been completely uprooted from the ground. 
Destroyed monuments usually can not be tied to previous surveys because of the high 
accuracy required for crustal motion research. If the suspected deformation is 
sig n i fican tly large, however, some information may be recovered if the monument (or a 
substitute) is reset in approximately the same location. Site inspection at MOUN clearly 
showed the position of the old monument and in early April, 1090, a rebar rod acting as 
a t em porary benchmark was set at approximately the same location. We estimate the 
temporary mark was set within 0.15 m of the previous monument. Because of its 
proximity to the 1987 rupture zone, the calculated seismic displacement at MOUN is 
about 0.5 m. Therefore, reoccupation of this site should retrieve a tectonic component 
larger than the expected uncertainty. 

The 1990 survey was conducted to establish the displacement of MOUN relative to 
its 1986 position. This "mini-campaign" included stationary GPS receivers at only three 
sites (Table 1). It was performed interactively with kinematic GPS transects along the 
southern San Andreas fault. Data were collected for only one night; a total of 9 satellites 
were visible during the scheduled 7.0 hour experiment. 

GPS Processing 

The 1986 GPS observations were processed with the GPS22 software developed at 
the National Geodetic Survey. Satellite orbital information was provided by the NSWC 
(Naval Surface Weapons Center). A tropospheric delay parameter was solved at each 
station, constrained by surface meteorological measurements. Ambiguities were fixed to 
the nearest integer. Each of the 20 days of observation was processed separately, and 
the daily solutions were combined to form one set of station coordinates with the 
geodetic adjustment program DYNAP (DYNamic Adjustment Program) [Drew and Snay, 
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1080]. All coordinates were computed in the WGS-84 reference frame [Defense Mapping 

& 

Agency, 1087]. Due to the relatively poor data quality, the accuracy of these 
measurements is believed to be on the order of 1 ppm (parts per million). 

Data from both 1088 campaigns were processed with the Bernese GPS analysis 
software (version 3.0), from the University of Bern in Switzerland. (Preliminary 
processing such as data translation and cycle slip fixing were performed with an earlier 
version of the software.) The Bernese GPS package allows 3-dimensional station 
coordinates to be determined from the integration of carrier phase, pseudo-range, and 
orbital data ]e.g., Gurtner et al . , 1085; Beutler et al., 1085; Rocken, 1088]. In addition to 
measuring station coordinates, the double difference estimation algorithm can solve for 
adjustments to six Keplerian elements and two radiation pressure coefficients for each 
satellite, an atmospheric zenith delay parameter at each station, a clock error term at 
each site, ionospheric model coefficients, and cycle ambiguity terms (e.g., Rocken, 1088]. 

For each of the 1088 surveys, all data were combined into a simultaneous multi-day 
solution. Surface meteorologic data (temperature, pressure, relative humidity) were used 
to constrain a Saasmotonian atmospheric model, and independent tropospheric zenith 
delay parameters were estimated at each station. We experimented with fixing 
ambiguities but found mixed results, therefore, ambiguities were left unresolved in the 
final analysis. 

In addition to station coordinates, satellite orbital parameters were estimated for 
both solutions. GPS observations from the CIGNET tracking sites were used to 
constrain the orbits. The Bernese software is able to combine ephemeris information 
from several days into a single multi-day arc. It has been suggested that multi-day 
satellite arcs significantly improve GPS precision [e.g., Lichten, 1887]. Typically, 6 
Keplerian elements and 2 radiation pressure coefficients are estimated for each satellite. 


11 


Davis ct al. [1989] utilized 4-day arcs (72 hours plus the length of the daily observation 
session) to analyze GPS data from North America. A 5-day satellite arc was used to 
process GPS data off the southern California coast near Santa Barbara [Larsen et al ., 
1991]. 

There comes a point, however, when the force model describing the satellite orbits 
begins to break down. It is unlikely that the 10-day interval of the February /March 
1988 survey (9 observation days plus one day of no measurements) could be modeled 
from a 10-day satellite arc. For purposes of orbit determination, therefore, this 
campaign was divided into multiple 3-day arc segments (Julian days 55-57, 58-60, and 
62-64), each defined by its own initial conditions. The data were still processed 
simultaneously, except that independent satellite parameters were estimated for each 3- 
day arc. For the entire campaign, 18 Keplerian orbital elements and six radiation 
pressure terms were estimated for each satellite. A similar technique was used for the 
March/April survey. The satellite orbits for this campaign were defined by two 3-day 
arcs (Julian days 88-90 and 91-93). 

Orbits were improved by holding the coordinates of at least 3 CIGNET stations 
fixed, to values well determined from VLBI and satellite laser ranging [Murray and King, 
MIT interoffice memorandum, 1988]. The GPS phase observables from these fiducial 
sites constrain the orbits, which in turn improve solution accuracy of the which in turn 
improve solution accuracy of the unfixed stations. There was always at least one day in 
each 3-day arc interval where observations from 3 fiducial sites were available [Larsen, 
1991]. 

After each 1988 campaign was processed separately, the Cartesian coordinate 
differences from the two solutions were adjusted by least squares to obtain one set of 
coordinates for that year. 
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The 1990 "mini-campaign" was processed with the Trimvec software, made available 
from the receiver manufacturer (Trimble Navigation, Sunnyvale, California). Recall that 
the purpose of this survey was to establish the displacement at station MOUN, presumed 
large because of its proximity to the 1987 seismic rupture. Since this mark had been 
reset between the 1986 and 1990 occupations, the error in the displacement estimate will 
be fairly large (~ 15 cm). Therefore, high accuracy from a GPS perspective was not 
required. The orbits were given by the broadcast ephemeris, and surface meteorological 
data were used to constrain a tropospheric delay model. 

Station Displacements 1986-1988 

GPS station displacements for the interval 1986 to 1988 are shown in Figure 4. All 
movements are made relative to station OCTI. Only the horizontal components are 
shown. The method in which errors are formulated and utilized is disussed below. 
Generally, the observed displacements can be decomposed into 3 components: 1) seismic 
deformation due to the Superstition Hills earthquake sequence, 2) secular deformation 
due to the Pacific-North American relative plate motion, and 3) unusually large 
measurement error, most apparent in the east-trending direction. 

The GPS displacement vectors suggest considerable deformation between the 1986 
and 1988 campaigns, of which a significant fraction is attributed to the 1987 Superstition 
Hills earthqu ake sequence. Stations nearest the seismic rupture zone (KANE and L589) 
show movements on the order of 40 cm. In fact, the 13 km KANE-L589 baseline was 
shortened by 70 cm. The orientations of the displacements are consistent with the 
conjugate fault pattern indicated by the mapped surface ofisets (i.e., right-lateral rupture 
on the Superstition Hills fault and left-lateral rupture along the Elmore Ranch fault). 
Other stations near the active fault system appear to have been affected by the 1987 


event. 


- 13 - 


Xhere is an additional component of displacement not readily explained by the 
seismic deformation. Stations east of the Imperial fault tend to be moving south or 
southeast relative to sites on the other side of the valley. This secular displacement is 
attributed to the relative motion between the North American and Pacific plates. 
Unfortunately, it is difficult to ascertain the magnitude of this deformation since station 
coverage west of the Imperial fault is somewhat lacking and many of these sites have 
rather large seismically induced displacements (recall that all GPS vectors are relative to 
OCTI). However, the deformation does appear to be considerable, which is not too 
surprising considering conventional geodesy indicates a significant fraction of the relative 
plate motion is occurring across the Imperial Valley. 

Also evident in Figure 4 are unusual movements which do not appear to be 
tectonically related. Most notable are the southwest-trending vectors (as opposed to 
southeast) for those sites near the border east of the Imperial fault. It appears as if the 
entire network has undergone a systematic clockwise rotation. We have investigated this 
possibility by assuming the network could be rotated (and translated) in terms of an 
outer coordinate solution by minimizing the displacement component perpendicular to 
the structural axis of the valley (N40*W) [e.g., Prescott, 1981]. Stations KANE and 
L589 were not included in the solution. However, the applied adjustment did not correct 
for the anomalous displacements, and in fact, made the apparent deformation less 
uniform. Although these unusual movements can not be attributed to a simple 
coordinate rotation, they can be explained by large east-trending systematic errors in the 
1986 data. This is consistent with the longitudinal orientation of the computed error 
ellipses (see below), and suggests the north-trending displacement components may be a 
more reliable indicator of tectonic deformation. 


Station Displacements 1086-1090 
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The 1986-1990 displacement of MOUN relative to BLAC is shown in Figure 5 

t. 

(dashed arrow). Although the uncertainty is large because MOUN was reset between 
surveys, the GPS data indicate significant movement attributed to the Superstition Hills 
events (recall MOUN is less than 1 km from the 1987 surface rupture). However, based 
on conventional geodesy as well as the 1986-1988 GPS movements, a non-seismic 
displacement component is suspected in the measurement. We attempt to remove this 
component by estimating the MOUN-BLAC secular displacement based on EDM 
observations. This is discussed in more detail below. We remove 2 years of the 
accumulated MOUN-BLAC secular movement, and then compute the displacement 
relative to station OCTI. Although the measurement still contains 2 years of non- 
seismic deformation, it places the displacement into the same reference frame as the 
1985-1988 movements. The adjusted MOUN displacement is shown in Figure 5 (solid 
arrow). 

Station Displacements 1988-1988 (February-April) 

Seven GPS sites were occupied during both 1988 campaigns; calculated site 
displacements for this 1 month interval are shown in Figure 6. An adjustment (simple 
translation) was applied to all movements so that the sum of the vector displacements is 
zero. The magnitude of the apparent movements range from 0.9 to 2.9 cm, averaging 
1.8 cm. It is interesting that the 2.9 cm displacement at KANE is to the west-southwest. 
This would be expected if left-lateral afterslip occurred along the Elmore Ranch fault. 
However, while postseismic oflsets were significant along the Superstition Hills fault, 
almost all Elmore Ranch activity ceased after the initiation of the 2nd main event 
[ Williams and Magistrate, 1989; Magistrate, 1989; Hudnut et at., 1989a]. The observed 
displacements in Figure 6 probably represent measurement error and not real 
deformation. Note that the largest vector components are oriented along the longitudinal 


direction. Station 0217 also exhibits a fairly large apparent east-trending movement, 
although this is almost certainly not tectonically related. 

Even if the displacements shown in Figure 6 are entirely measurement error, they do 
illustrate two points: 1) GPS accuracy is easily sufficient to monitor tectonic motions, 
and 2) the anomalous displacements for the 1086-1088 interval (Figure 4) are probably 
due to the quality of the 1086 data. Typical crustal deformation rates across major 
tectonic structures in southern California are between 1 and 5 cm/yr. At 1 to 2 cm GPS 
accuracy, such deformation might be resolvable in as short as 1 year. Although the 
short baseline (~ 50 km) horizontal precisions computed from GPS repeatability studies 
are generally at the sub-centimeter level [e.g., Dong and Bock, 1080], these tests usually 
involve multiple occupations of the same network over a consecutive 4 to 5 day interval. 
The 7 stations shown in Figure 6 were not all observed simultaneously [Larsen, 1001]. 
Therefore, the repeatability is somewhat degraded since errors tend to propagate through 
the solutions. However, the relatively good consistency suggested with the 1088 results 
suggest that some of the unusual movements observed with the 1086-1088 displacement 
vectors are probably due to the quality of the 1086 data. 

GPS Errors 

Formal estimates of GPS uncertainty almost always underestimate variances derived 
from repeatability studies. We attempt to determine more realistic and illustrative 
errors by multiplying the formal covariance matrix calculated with each GPS solution by 
an estimated variance factor, which scales as the average baseline length. We have 
chosen a variance factor for the 1086 data so that the average baseline error scales as 1 
ppm. Since the 1088 data were processed robustly utilizing orbit improvement 
techniques (not available with the 1086 data), we assume these errors are negligible. 
Although this method is somewhat ad hoc, it does illustrate an important fundamental. 
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Notably, the east-trending error is about 4 times larger than the north-trending error. 
This is attributed to the north-south ground-track of the satellite orbits. 

Because station MOUN was reset between the 1986 and 1990 surveys, the largest 
error source is due to the difficulty in establishing the new mark at the previous location. 
From field inspection, this uncertainty is estimated to be 15 cm (in all directions) and is 
additive to the GPS measurement error. 

5. Modeling 
Theory 

Simple dislocation theory is often used to model seismically induced geodetic 
deformation. The earth is considered a homogeneous isotropic elastic half-space with no 
stress applied to the free surface. The displacement field u k for a dislocation E in the 
medium is given by 

a) 

where Vuj is the discontinuity, w* are the displacement Green’s functions due to a set of 
strain nuclei, and Vj are the direction cosines of the normal to the surface element dE 
[Steketcc, 1958; Chinnery, 1961]. Analytical solutions to this integral are rather complex, 
but have been simplified for special cases of dislocation or fault geometry [e.g., Chinnery, 
1961; Savage and Hastie , 1966; Mansinha and Smylie, 1967]. General expressions of the 
displacement field fpr rectangular strike and dip-slip faults of arbitrary Inclination have 
been computed by Mansinha and Smylie [1971] and Okada [1985]. Arbitrary slip 
directions can be designed by the superposition of strike and dip-slip dislocations. 

The strain/stress within a medium is computed by differentiating the displacement 
field. For the displacement u k , where ti is a function of the geometrical coordinates x if 
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the strain tensor components are given by 




chi; 5U; 

+ 1 


( 2 ) 


,J 2 l&Cj dx> j 

In an isotropic medium the stress tensor er is given by 

o-j-X^j+2/xeij 

s 

where 0 is the dilatation (0— ]T) c ii)- We ass'™* the medium is a Poisson solid with 

i-i 

2.8X10 11 dyne-cm -2 . 


( 3 ) 


A rectangular dislocation within an elastic half-space will create a spatially 
dependent stress tensor <r throughout the volume. The force acting at a point along an 
arbitrarily oriented plane in the medium is computed by multiplying a by the outward 
normal vector to the plane (ft n ). That is, the traction vector ¥ is given by 


¥«of$ n (4) 

If we assume the plane is coincident with a fault, then the forces generated on this 

secondary structure due to the initial dislocation are determined by calculating the 

traction vectors at selected points along the fault. The normal (cr n ), strike-slip (<r 8 ), and 

dip-slip (cr d ) stresses on the fault plane are computed by 


<r n «¥ft n 

(5a) 

<7 8 «T-fi 8 

(5b) 

a d -¥-fl d 

(5c) 


where fl n , K 8 , and ft d are the normal vectors perpendicular, along strike, and along dip 
to the fault plane. Analytic solutions for the stress and strain fields due to a dislocation 
within a medium are given by Iwasaki and Sato [1070] and Altxmnt [1074]. 


Inverse Methods 


Inversion of seismically generated geodetic displacements can yield fault-rupture 
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parameters, such as the slip distribution along a fault plane [e.g., Ward and Barrientos, 
1986; Harris and Segall, 1987; Segall and Harris, 1987; Snay , 1989]. We employ a method 
similar to that outlined in Segall and Harris [1987]. Singular Value Decomposition (SVD) 
[e.g., Lanczos, 1961; Jackson, 1972; Menke, 1984] and elastic dislocation theory are used 
to invert the Imperial Valley GPS measurements for seismic slip along the Superstition 
. Hills and Elmore Ranch faults. 

The relationship between surface deformation and slip along a rectangular dislocation 
is defined by Equation (l). The rupture plane is modeled as a set of non-overlapping 
rectangular dislocations. That is, the fault plane is partitioned into multiple sub- 
elements or sub-faults. The slip distribution along the seismically active fault is given as 

* 

the discreate approximation of slip along each sub-element. The normal equations which 
govern surface displacement resulting from such slip is given by 

A g m f -* d g (6) 

where the superscripts g and / refer to geodetic observation and fault slip, respectively. 

Each row of A g is determined from (1), and is a function of sub-fault geometry and 

geodetic position. The slip distribution m f is defined by m f — [m 1 ,m 2 , • • • ,m n ] T , where 

nij is the slip along the rth sub-fault. The data vector d 8 contains the geodetic 

observables. Each GPS station displacement will add 3 rows to A 8 and 3 elements to d 8 , 

corresponding to the vertical and two horizontal components. In practice, we may 

choose to ignore the less accurate vertical observation. This is especially true in strike- 

slip environments where horizontal displacement predominates. 

Surface rupture is easily included into (6) by considering measurements of surface 
displacement as geodetic observation. The ofisets are modeled as a priori slip 
information on the surface intersecting sub-faults. Equation (6) becomes 
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•where d® are the discreate approximations of surface slip along the fault trace, and 
A*-l if sub-fault element j corresponds to surface slip ofiset », otherwise Ay — 0. 


The GPS displacements shown in Figures 4 and 5 are not connected to an external 
reference but are defined relative to station OCTI. (Some geodetic measurements such as 
line-length change recorded by EDM observations are Independent of an absolute 
reference.) The displacement at OCTI is assumed to be 0, which in fact may be true from 
a seismic standpoint since this site is far from the Superstition Hills rupture zone. 
However, any attempt to use the GPS displacements as a criteria for evaluating the effect 
of the earthquake sequence will be distorted by measurement error at OCTI and/or non- 
seismic deformation between this site and the other stations. This ambiguity is largely 
circumvented if displacement-offset terms are estimated in addition to the fault slip 
parameters. Equation (7) is then rewritten 


Am-A 


m‘ 


nr 


«d 


( 8 ) 


where m® is the tth non-seismic component uniformly added to all station displacements 
(i.e., north, east, vertical). 


The Singular Value Decomposition of A is given by 

A-UXV 1 (9) 

where U is a matrix of eigenvectors spanning the data space, V is a matrix of 

eigenvectors spanning the parameter space, and X is a diagonal matrix of singular values. 

Without loss of generality, this is written 

A-UpXpVT (10) 

■where p refers to the non-zero singular values. If the normal equations of (8) are 

normalized to have unit variance [e.g., Scgall and Harris , 1987], the generalized inverse of 
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(8) and (10) is given by 

A^-VpX^Up (11) 

[Lanczos, 1961; Mcnkc, 1984]. In practice it is often necessary to restrict the volume of 

the parameter space by considering only the k largest singular values, setting all others 

to 0. 

The generalized solution to (8) for the k largest singular values is given by 

m “A^d+VoOfo (12) 

where V 0 are eigenvectors spanning the null space of the model and is a vector of 

arbitrary coefficients. The volume of the model space not constrained by observation is 

defined by V 0 c^. This term is not influenced by the geodetic data and is thus arbitrary. 

■ Often it is the minimum-length solution m ■»A 1 f 1 d which is of interest (the coefficients of 
c*b are 0). However, some other solution criteria can be satisfied by carefully designating 
the coefficients of a^. 

For high resolution fault models where the rupture plane is partitioned into 
numerous sub-faults, it is necessary to apply some type of smoothing constraint over the 
dislocation surface to prevent the slip distribution from taking on an oscillatory pattern. 
ScgaU and Harris [1987] showed that the "roughness" of fault slip could be minimized by 
considering smoothness as an a priori constraint utilized from the model null space 
through the coefficients of ocq (Equation 12). They used a smoothing matrix T, with 
coefficients determined from the discreate approximation of the Laplacian operator 
^ — + where x and y are the fault distances along strike and dip, 

respectively. The boundary conditions around the lower and lateral edges of the 
dislocation are assumed to be null slip, so that the applied smoothing operator causes the 
calculated fault offset to tend toward zero along these boundarys. Because the 
Superstition Hills and Elmore Ranch faults ruptured the surface, the upper boundary is 
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considered an unconstrained dislocation. The estimated fault slip is given by 

in - |I - V 0 (V 0 t T t TV 0 )- 1 'Vo^TjAk-M (13) 

(Equation 13, [Segall and Harris, 1087]). A similar formulation for utilizing fault 

smoothness over the model null space is given by Harris and Segall [1987]; an alternate 

method considering fault smoothness as quasi-data is provided by Snay [1089]. 

For over-constrained solutions, where there are more independent data than 
parameters, if k — p then SVD is equivalent to simple least-squares. This is 
advantageous since the solution provided by (12) can be utilized for either uniform 
dislocations or for detailed parameterizations where the fault plane is partitioned into 
multiple sub-elements. 

Simple dislocation theory has the advantage that the displacement and stress/strain 
fields for simple fault ruptures can be computed almost instantaneously. The inverse 
problem of using geodetic data to calculate the slip distribution along a fault plane is 
also straightforward. However, we have assumed the earth can be modeled as a 
homogeneous half-space. Crustal layering or inhomogeneities in the earth can introduce 
non-existent structure into half-space models [ Savage , 1087]. While low-resolution 
schemes such as the average slip over the fault plane will not be seriously affected, 
attempts to resolve detailed properties may be badly contaminated by artifacts of earth 
structure. 

Seismic Displacement 

We model rupture along the Superstition Hills and Elmore Ranch faults as strike-slip 
dislocations along vertical planes extending from the surface to 10 km depth. Each 
dislocation approximately coincides with the mapped surface rupture and/or aftershock 
distribution. The geometrical parameters for the modeled faults are listed in Table 2. 
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Initially, the Superstition Hills and Elmore Ranch faults are modeled as uniform 
dislocations, with no slip variation allowed on the rupture planes. The dislocations are 
then partitioned into multiple sub-elements and the slip distribution along the two faults 
is calculated from the discrete offset for each sub-fault. 

For an initial estimate of slip along the two faults, we use only the displacements at 
KANE, MOUN, and L589; the three GPS sites nearest the seismic rupture zone (Figure 
7). Because the observed movements at these stations are relative to OCTI, we also solve 
for uniform north and east-trending offsets in the displacements. This will remove any 
systematic distortion due to measurement error at OCTI. Both faults are regarded as 
uniform dislocations. The solution to this simple model is given in Table 3 (Model 1) 
and the observed vs. calculated displacements are shown in Figure 7. Generally, there is 
good agreement between model and observation. This is not altogether unexpected since 
we are solving for 4 parameters with 6 data. The large residual at MOUN is presumably 
due to the uncertainty in relocating the 1990 reset monument at the 1986 position (recall 
the estimated error is 15 cm in each direction). The surface deformation pattern 
computed from Model I is fairly extensive; the 4 cm horizontal displacement contour is 
shown in Figure 7. 

The north and east-trending displacement components at each Imperial Valley site, 
plotted as a function of distance from a N40 * W trending line through OCTI are shown 
in Figure 8. Shown are the non-selsmic movements; that is, the seismic component 
computed from Model 1 is subtracted from the observed displacements (Figure 4). 
Stations where the seismic correction is greater than 4 cm (in the component plotted) are 
shown as open circles, other sites as filled circles. The displacements represent the non- 
seismic deformation perpendicular to the plate motion. 

A fairly well defined pattern is observed in the north-trending components (Figure 8). 
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Stations to the northeast have moved south (or southeast) about 8.1 cm relative to sites 
on the other side of the valley. Stations which display the largest deviation are for the 
most part those sites where the applied seismic correction Is greater than 4 cm (open 
circles). This may indicate fault complexity not accounted for by the simple dislocation 
used to model the 1987 earthquake. The east-trending displacements, however, show 
large data scatter; no distinguishable pattern is readily visible across the valley. The 
scatter is not a function of the magnitude of the applied seismic correction, so 
presumably it is due to the large measurement errors inherent in the longitudinal 
direction. 

We assume the 8.1 cm north-trending offset is attributed to plate-boundary 
deformation between the North American and Pacific plates (Figure 8). Taking into 
account the orientation of the suspected deformation (N40*W), as well as the time 
interval between the 1986 and 1988 surveys (1.79 years), the north-trending movements 
suggest 5.9 cm/yr displacement across the valley. This is significantly larger than the 
3.45 - 4.3 cm/yr rates obtained from conventional surveys. Although accelerated 
deformation between the GPS campaigns can not be ruled out, there is relatively poor 
station coverage in the southwest portion of the valley so it is difficult to precisely 
estimate the valley crossing displacement. This is especially true considering most of the 
southwestern sites suffered large seismic displacements during 1987. If Model 1 (Table 3) 
does not accurately reflect the rupture process during the Superstition Hills earthquake 
sequence, the unmodeled effects will influence the non-seismic estimate. 

The seismic rupture estimated in Model 1 depends heavily on shallow slip since the 
three stations used are within close proximity to the dislocation planes. While field 
mapping shows that rupture extended to the surface, slip at depth may go undetected. 
Therefore, it is necessary to examine the displacement at stations away from the fault to 
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ascertaln the extent of faulting with depth. 

To incorporate more GPS data into the fault-plane inversion, it is necessary to 
remove the non-seismic deformation from the displacement field. Perhaps the best 
example of secular deformation across the Imperial Valley is provided by U.S.G.S. 
trilateration measurements between 1972 and 1987 [e.g., Prescott ct al., 1987a; Prescott 
et aL, 1987b]. Computed station velocities for the Salton Trough EDM network are 
roughly parallel to the direction of plate motion (N40*W) (Figure 9), although to some 
degree the geodetic orientation is dictated by the outer coordinate solution imposed to 
transform EDM line-length changes into station displacements [Prescott, 1981]. The 
total differential velocity across the network is 3.45 cm/yr and is accommodated in a 50 
km wide zone [Prescott et al., 1987b]. Triangulation measurements suggest a larger rate 
i between 1941 and 1986 (4.3 cm/yr), but these observations are not as accurate as the 
EDM measurements [5nay and Drew, 1988]. However, along the Imperial fault where 
* EDM sites are sparse, the triangulation data suggest concentrated deformation in a 
r narrow 20 km wide zone. The conventional geodetic data are modeled using the 
following empirical approach. The differential velocity across the valley is assumed to be 
t 3.45 cm/yr. Running along the axis of the valley is a transition zone(s) (Figure 9), where 
the strain gradient is defined by simple shear with the displacements oriented N40 * W. 
The transition zone north of the Imperial fault Is 50 km wide; to the south it is 20 km 
wide. The modeled station movements fit the observed EDM displacements extremely 
well (Figure 9). 

The secular deformation derived from the conventional measurements is subtracted 
from the observed GPS displacements (Figures 4 and 5), leaving the seismic component 
(and measurement error). For Model 2 (Table 3), the uniform slip along the Superstition 
Hills and Elmore Ranch faults is recomputed using all GPS data (without the secular 
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deformation). The residuals (observed minus calculated) for this model are shown in 
Figure 10. The largest residuals trend in the longitudinal direction and are especially 
noticeable for those sites in the southeast. This reconfirms our speculation for large 
east-trending errors. However, the residual vectors at the three stations nearest the 
seismogenic zone are unusually large. The large vector at MOUN is easily explained since 
this station was reset between surveys. The residuals at L589 and KANE, however, are 
significantly larger than the average residual computed at the other stations. Field 
inspection at both localities revealed no evidence that either monument had been 
disturbed. Because both sites are located near the seismic rupture zone, this suggests 
additional slip not accounted for by the simple dislocation parameters used for Model 2. 

Although the estimated slip is marginally different between Models 1 and 2 (Table 3), 
this is not significant considering the estimated uncertainties. In fact, because the near- 
field (Model 1) and far-field (Model 2) solutions are similar, this suggests that to first 
order there is not significant slip dependence with depth (within a factor of 2 or 3). It is 
alsr> noteworthy that the uncertainties improve by only ~ 50 % with the additional data 
supplied with Model 2. This illustrates the necessity of measurements near the seismic 
rupture. For simple fault models, where uniform slip is constrained to a single 
dislocation plane (or two planes in the case of the Superstition Hills sequence), it is as 
important to have at least minimal station coverage within a few kilometers of the 
seismogenic zone, as it is to have many sites located away from the fault(s). 

'** Seismic Slip Distribution 

Uniform dislocations for the Superstition Hills and Elmore Ranch faults were utilized 
for Models 1 and 2. To estimate the seismic slip distribution, it is necessary to partition 
the rupture planes into multiple regions or sub-faults (Model 3). The divisions must be 
sufficiently dense as to provide reasonable slip resolution. We choose 10 sub-fault 
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elements in the horizontal and 5 in the vertical, so that each fault is partitioned into 50 
sub-regions. The dimensions of each dislocation element is 2.0 km in width (vertical) 
and about 2.5 km in length (2.4 km for Superstition Hills and 2.6 km for the Elmore 
Ranch fault). The slip distribution along the two rupture planes is constrained to be 
sufficiently smooth (Equation 13). The GPS displacements are adjusted according to the 
estimated north and east-trending offsets estimated for Model 2. 

In addition to the GPS data, a priori surface-slip information is added to the 
solution. Surface slip along the Superstition Hills fault extends (nearly) the entire length 
of the modeled fault plane (Figure 11) [Williams and Magistrals, 1989]. The surface 
rupture has been incremented into 2.4 km segments corresponding to the horizontal 
dimension of each sub-fault. The average slip over each segment is assigned as an a 
priori slip estimate for the surface fault element to which it corresponds. Surface rupture 
along the Elmore Ranch fault is confined to the southwestern segment (Figure 11). 
Recall for this event that the mapped surface breaks occurred along several nearly 
. parallel strands (Figure 2) [. Hudnut et al., 1989a]. We take the cumulative surface offset 
for all strands averaged over 2.6 km segments, and use this as a priori slip for the 
surface sub-fault elements. Where no rupture is mapped (to the northeast), the surface 
t intersecting fault partitions are assigned 0 slip. The a priori uncertainty for each 
surface-slip estimate is assumed to be 10 cm. 

The number of independent parameters estimated through singular value 
decomposition depends on the number of singular values k utilized In Equation (12). A 
trade-off ousts between solution variance and resolution [e.g., Mtnkt, 1984]. While large 
k produces highly resolved models, this is at the expense of increasing solution 
uncertainty. Correspondingly, small k yields low variance solutions but does not provide 
detailed resolution. A total of 100 sub-fault elements are use for Model 3 (50 for each 
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fault). If k — 100, then slip along each sub-fault will be determined uniquely. Because of 
limited geodetic coverage, however, it is practical to consider only the first few 
eigenvectors of the parameter space defined by the geodetic observations. Therefore, each 
solved parameter Is a function of some combined slip over multiple sub-fault elements. 
This is a fundamental property of singular value decomposition when used to solve 
under-determined or poorly-determined problems [e.g., Jackson, 1072]. It is necessary to 
determine the k which maximizes the resolution without allowing the solution to become 
too oscillatory or unstable. 

The geodetic moment, solution instability, and model residual calculated for different 
values of k are shown in Figure 12. The moment is a function of the average slip along 
the fault planes, while solution instability is determined from the standard deviation of 
slip for each sub-fault element. An instability of 0 (stable) indicates uniform slip along 
the fault planes, while high values indicate an oscillatory or unstable solution. The RMS 
indicates the agreement between model and observation and is calculated by 
RMS — $](oi — cfi/a? where Oj is the observed, q is the calculated, and a- is the 
uncertainty assigned to the ith observation. 

When surface slip is incorporated into the solution (Figure 12a), the first 20 singular 
values are well constrained by the measured offset along the fault, and are influenced 
little by the geodetic observations. The surface measurements reflect the dislocation on 
the uppermost fault elements, with little depth resolution. This is illustrated by the 
solution for k — 20. The calculated moments are significantly less than that for Model 2, 
pres um ably because the surface offset is not representative of the larger displacement 
along the rest of the fault plane(s). Consequently, in order to estimate the slip 
distribution with depth it is necessary to consider solutions where k > 20. The solution 
becomes increasingly oscillatory after k * 30 as is indicated by the instability. The RMS 
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is significantly reduced beyond k — 20, but only improves marginally -with increasing k. 
The solution fit for k > 20 is slightly better than that for Model 2. 

Also evident in Figure 12a are the large moment estimates for the Elmore Ranch 
fault, ■which are almost equal to the computed moments for the Superstition Hills fault. 
This is unexpected considering the latter event yielded a significantly greater surface 
wave magnitude, as well as a larger moment estimated from Model 2. It is difficult to 
distinguish slip between the two faults using SVD, as is illustrated by Figure 13. The 
displacement pattern for a 130 cm rupture on the northwest-trending right-lateral 
Superstition Hills fault is compared with a 30 cm rupture on the northeast-trending left- 
lateral Elmore Ranch fault. The deformation pattern between the two dislocations is 
almost identical (although different in magnitude). As a result, the larger slip estimate 
along the Superstition Hills fault is being mapped onto the Elmore Ranch fault plane 
producing the higher moment. This is especially true since only the first few model-space 
eigenvectors are used, as is necessary because Model 3 is very underdetermined (more 
model parameters than data). Therefore, only linear combinations of model parameters 
are uniquely defined. The similar moments indicated by Figure 12a suggest slip between 
the two faults is strongly correlated. 

The estimated seismic slip distribution along the Superstition Hills and Elmore 
Ranch faults for k — 23 and k — 27 are shown in Figure 14a. We refer to these solutions 
as Models 3a and 3b, respectively (Table 3). Although both faults are partitioned into 
50 elements, the fault-rupture is not as resolved as the contours suggest since only the 
first few (non-surface) model-space eigenvectors are independently solved. The solution 
suggests fairly uniform rupture along both fault planes for k * 23. The dislocation may 
be slightly concentrated to the southwest along the Elmore Ranch fault. The apparent 
"bullseye" pattern is due to the smoothness constraint, which requires zero slip along the 
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lateral edges and lower boundaries (the upper boundaries are constrained by the surface 
slip). There is little difference in the slip distribution for k-21 through k-26. 
However, there is a noticeable change in the dislocation pattern starting with k — 27. 
While slip along the Elmore Ranch fault still appears fairly uniform, displacement along 
the Superstition Hills fault is concentrated to the northwest and to the southeast. This 
change is significant and is caused by the GPS displacement at one station. Recall the 
large residual for L589 in Model 2 (Figure 10). This discrepancy is nearly eliminated 
beginning with k —27. Therefore, in order to satisfy the observed displacement at L589, 
it is necessary to concentrate rupture at each end of the Superstition Hills fault. This 
analysis assumes the observed GPS displacement at L589 is seismically generated, and 
not contaminated by an unusually large measurement error. The slip deficit near the 
center of the fault roughly corresponds to the drop in slip measured at the surface 
(Figure 11). 

Independent solutions are made without constraining the upper sub-fault elements by 
surface offset measurements (Figure 12b). The unconstrained moments are generally 
smaller than when surface slip is used. This is because the surface measurements are less 
than the average slip estimate along the fault plane. Without surface constraint, the 
geodetic data are satisfied to a greater degree by slip near the surface; otherwise, it is 
necessary to compensate the small shallow offsets by increased slip at greater depths. 
The seismic slip distribution estimated without surface constraint along the Superstition 
Hills and Elmore Ranch faults for k — 3 and k— 7 are shown in Figure 14b. We refer to 
these solutions as Models 3c and 3d, respectively (Table 3). For these unconstrained 
solutions, the slip distribution is confined to shallower depths. Nevertheless, the solution 
for k— 7 yields concentrated slip at each end of the Superstition Hills fault. This 
emphasizes that the calculated slip discrepancy near the faults midsection is not simply 
an artifact of incorporating measurements of surface offset into the model. 
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6. Geophysical Implications 
Conjugate Faulting 

The most prominent feature of the Superstition Hills earthquake sequence is the 
conjugate relationship exhibited by near-simultaneous ruptures along right-lateral 
northwest and left-lateral northeast-trending faults. In the context of the Imperial 
Valley, the northeast-trending structures are termed "cross-faults" [e.g., Hudnut it al., 
1080a]. Conjugate and cross-fault seismicity seems to be a fairly typical phenomenon for 
this region (Figure 15), and may dictate the strain transfer mechanism between faults. 
The 1081 Westmorland earthquake (Ml 4.1) is a prime example of cross-fault tectonics. 
The mainshock and aftershock sequence is clearly mapped onto a northeast-trending 
lineament. Other examples are associated with the Imperial fault. The largest 
aftershock (Ml 5.8) following the 1070 Imperial Valley earthquake (Ml 6.6) was located 
near the town of Brawley [Johnson and Hutton, 1082]. The focal mechanism and 
following seismicity suggested left-lateral slip along a vertical northeast-trending fault. 
Riilinger and Larstn [1086] found that rupture along an identical conjugate structure 
successfully modeled geodetic observations within the Brawley Seismic Zone. A large (Ml 
5.5) aftershock was also recorded near Brawley following the 1040 earthquake [Neumann, 
1042]. Due to the sparsity of seismic data, neither the mechanism nor location are 
precisely determined, although we speculate this event occurred along the same 
northeast-trending feature as the large 1070 aftershock. Of historical interest are 
Imperial Valley earthquake pairs during 1015 (Ml 6.3, Ml 6.3) and 1027 (Ml 5.8, Ml 5.5) 
[Beat, 1015; Toppozada it al., 1078]. In each case, the 2nd shock followed the first by 
about 1 horn, contrasting with the 12 hour interval between the 1087 events. It is not 
known which fault(s) ruptured during these earthquake sequences, but conjugate fault 
interaction is highly probable. 
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Rupture on the Superstition Hills fault was almost certainly triggered by the Elmore 
Ranch event (occurring 12 hours earlier) suggesting some mechanism of stress transfer 
between the two faults. Figure 16 shows the normal (<r„) and strike-shear (<r s ) stress 
components instantaneously applied to the Superstition Hills fault due to a 30 cm left- 
lateral Elmore Ranch dislocation (Model 2). Tension and right-lateral shear are 
considered positive, both tending to induce failure on the rupture plane. Also shown is 
the Coulomb failure stress (<r c ), here given by <r c —<r B +na n , where 0.75. Positive 
values indicate stress-loading leading toward shear failure. 

The stress regime necessary for left-lateral rupture along a northeast-trending fault is 
identical to that required for right-lateral failure along a northwest-trending fault. 
Hence, we can assume that the Superstition Hills rupture plane was at or near failure at 
the time of the Elmore Ranch event. The initial shock generated an increase in the 
Coulomb failure potential along the Superstition Hills fault (Figure 16), possibly 
advancing it past its failure threshold. This is seen mostly as a combination of reduced 
compressive stress (earthquake inducing) countered by left-lateral shear (earthquake 
inhibiting). The increase is maximized along the northwest boundary of the rupture 
plane, near the nucleation point of the second event. Presumably rupture began where 
the applied stress was greatest and then propagated to the southeast. Northwestward 
rupture is prohibited because the increase in compressive forces tends to inhibit shear 
failure along this segment of the plane. The magnitude of the Coulomb stress increase 
near the Superstition Hills epicentral zone (~ 10 bars) is comparable to typical 
earthquake stress drops. 

The one to several hour delay recorded between events during observed and suspected 
conjugate episodes in the Imperial Valley is significant from an earthquake failure 
perspective. Shown in Figure 17 are potential scenario’s for earthquake ruptures 
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Involving conjugate-mainshock interaction, such as that observed for the Superstition 
Hills events. We assume faults fail by an undefined mechanism when they are at or 
above some critical stress level. The regional strain acting over several years brings a 
fault near this critical failure point. A stress increase is induced along part of the fault 
plane due to rupture on a conjugate structure (e.g., Figure 16), which may or may not be 
sufficient to push the stress state past its critical threshold. In the case of Earthquake 1 
(Figure 17a), the stress change caused by the conjugate event is not enough to induce 
failure. Some form of time dependent stress transfer onto the fault is activated and 
eventually the critical level is reached. Such a mechanism involving postseismic viscous 
creep along the Elmore Ranch fault has been suggested for the 1987 Superstition Hills 
sequence {Given and Stuart, 1988]. If this scenario is valid, we would equally expect 
failure modes such as that indicated by Earthquake 2. Here, the instantaneous stress 
applied to the fault from the conjugate event pushes the stress state past the critical 
level and rupture is immediate. In this case, failure along the two perpendicular fault 
planes will occur simultaneously. However, this behavior is not observed in the Imperial 
Valley. Conjugate episodes characteristically have been separated by one to several 
hours. This suggests that the critical stress level can be exceeded without immediate 
failure. Therefore, some time dependent mechanism must be active on the fault plane, as 
opposed to additional stress transfer through the crust. We loosely refer to this as 
"stress corrosion" (Figure 17b) [e.g., Das and Scholz, 1981]. In the case of Earthquakes 3 
and 4, it is suggested that the critical stress level must be exceeded for a period of one to 
several hours before failure occurs. Hudnut et al. [1989b] proposed fluid diffusion as an 
alternate mechanisnf, whereby the effective normal stress was reduced (made more 
positive) due to pore-fluid infiltration into the rupture plane, thus increasing the 
Coulomb failure stress. This process still involves action on the fault plane rather than 
stress tr ansf er. Regardless of cause, the temporal and geometric relationship exhibited 


- 83 - 


by the conjugate fault interaction is seemingly typical of Imperial Valley tectonics, and is 
thus an important factor for the potential prediction of large earthquakes and 
aftershocks. 


Moment and slip distribution 

The geodetic (GPS) source parameters for the Superstition Hills and Elmore Ranch 
earthquakes are listed in Table 3 and Figure 14. The seismic moment is best defined by 
Model 2, "while the slip distribution is best expressed by Models 3a and 3b. Model 1 is 
constrained with minimal station coverage and Models 3c and 3d do not include the 
additional information supplied from surface slip measurements. The GPS observations 
are directly proportional to the combined effect of the Elmore Ranch and Superstition 
Hills events, although we have attempted to resolve slip between each fault. The 
calculated parameters are a function of the coseismic offset, as well as 3 - 4 months of 
postseismic slip (plus 1.5 years of preseismic movement, if any). The average left-lateral 
Elmore Ranch dislocation is about 30 cm, with fairly uniform distribution along the 
fault plane. In the case of the Superstition Hills fault, the average right-lateral slip is 
about 130 cm, with concentrated deformation along the northwest and southeast sections 
of the fault. Because the GPS sampling frequency is so low (years), the calculated source 
parameters should contain the total coseismic moment, including several months of 
postseismic slip. 

The GPS moments are compared with estimates made through seismic and other 
geodetic studies (Table 4). Forward and inverse models using teleselsmlc and strong- 
motion recordings are used to constrain source parameters, as well as investigate 
complexities of the Superstition Hills rupture process [Dztcwonski et ai, 1080; Bent et ai, 
1080; Sipktn, 1080; Frankel and Wennerberg , 1080; Wald et al., 1000; Hwang et al., 1000]. 
The teleseismic moments agree fairly well with the GPS estimates, while the strong 
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ground motion data suggest significantly lower moments. The high-frequency strong- 
motion measurements are dominated by energy around 1 second and conceivably miss a 
sizable portion of the long-period energy release recorded with the teleseismic and GPS 
data. Hence, the near field seismic solutions may underestimate the total moment 
release. 

Geodetic measurements from Pinyon Flat observatoiy are used to constrain planer 
and curved dislocation models for the Superstition Hills and Elmore Ranch faults [Agnew 
and Wyatt, 1989]. The data are obtained from long-base strain and tilt-meters, as well 
as a borehole dilatometer. The best-fit planer models to all observations are significantly 
lower than those calculated with the GPS and teleseismic data (Table 4), although a 70 
% moment incr eas e for the Superstition Hills fault is obtained when the strainroeter data 
are excluded. The low moment estimate may be due to a number of factors [Agnew and 
Wyatt, 1989]: 1) measurement quality, particular with the strainmeter, 2) rheologic 
differences between Superstition Hills and Pinyon Flat, and 3) strainmeter-dilatometer 
sensitivity to the nodal deformation plane roughly on azimuth with the observatoiy. 

„ Geodolite observations from the Salton Trough EDM network were made in early 

v December (1987), several days after the two large earthquakes [Lisowski and Savage, 
1988]; the last previous occupation was in January, 1987. Simple dislocation models 
with 40 cm left-lateral slip along the Elmore Ranch fault and 120 cm right-lateral slip 
along the Superstition Hills fault best-fit the observations. The estimated moment for 
the Superstition Hills event is comparable to that obtained with the GPS displacements 
(Table 4). 

The discrepancies in Table 4 are largely attributed to the alternate methodologies, 
observations, and parameters used to constrain each model. However, for those 
calculations which include moment estimates for both the Superstition Hills and Elmore 
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Ranch events, the ratio between the two ruptures is fairly constant. This illustrates an 
internal consistency with each method. More importantly, it suggests that postseismic 
slip along the Superstition Hills fault is probably confined to the shallow segment of the 
rupture plane. Since seismic activity on the Elmore Ranch fault essentially ceased after 
the 2nd main event, if postseismic slip were occurring in mass along a large fraction of 
the Superstition Hills rupture plane, the GPS moment ratio would be considerably larger. 

While the epicenter and aftershock sequence for the Superstition Hills event were 
concentrated along the northwestern portion of the fault, seismic and surface slip data 
suggest significant moment release on the southern section of the Superstition Hills fault 
[Wald et al., 1989; Bent et al., 1989; Hwang et al. , 1990; Williams and Magistrals , 1989]. 
An exception is the strong ground motion study of Frankel and Wcnnerbcrg [1989], where 
slip is confined to the northwest. However, their low Superstition Hills moment (Table 
4) suggests this analysis may be strongly susceptible to the high-frequency content of the 
data, indicating rupture along the southeast segment was dictated primarily by low- 
frequency energy release. The GPS data also reveal dislocation along the southeastern 
segment of the fault, and further suggests a displacement null near the faults mid- 
section. This slip deficiency may be related to the surface oflset drop observed along the 
center of the fault (Figure 10). 

Deformation across the Imperial Valley 

The 1985-1988 GPS station displacements indicate a significant component of 
deformation across the Imperial Valley not associated with the 1987 Superstition Hills 
earthquake sequence (Figures 4 and 8). This motion is attributed to plate-boundary 
deformation due to the relative velocity between the Pacific and North American plates. 
From empirical evidence provided by Salton Trough EDM observations between 1972 
and 1987, these non-seismic movements were modeled as 3.45 cm/yr differential velocity 
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across the valley. However, after removing the seismic deformation predicted with a 
preliminary model (Model 1, Table 4) the calculated GPS displacements average 5.9 
cm/yr (Figure 8), significantly larger than that obtained with the trilateration 
measurements over the last two decades. This GPS velocity is also larger than the 4.7 
cm/yr plate motion predicted from global models (for Imperial Valley coordinates) 
[DeMcts et al., 1990]. The GPS rate may be an overestimate, especially considering the 
large measurement uncertainties and the influence from the 1987 earthquake sequence. It 
is posable, however, that the 1986-1988 GPS rate could represent accelerated 
deformation. In fact, GPS observations between 1988 and 1989 indicate 5.2 ± 0.9 cm/yr 
displacement across the Imperial Valley [Larsen and Reilinger, 1991]. Accelerated 
deformation is not without precedence. Triangulation observations suggest a rate of 6.2 
cm/yr between 1941 and 1954, although this is attributed to postseismic deformation 
following the 1940 Imperial Valley earthquake. Increased deformation following the 1979 
earthquake is not observed in the EDM observations et al -, 1986]. Additional 

GPS measurements are necessary in the Imperial Valley in order to ascertain the current 
deformation rate. 


7. Conclusions 

Station movements computed from 4 Imperial Valley GPS campaigns indicate large 
crustal displacements during the periods 1986-1988 and 1986-1990. Much of the 
deformation is attributed to the 1987 Superstition Hills earthquake sequence. Ten sites 
near the seismic rupture zone are displaced at least 10 cm, although the GPS 
observations contain large uncertainties. This Is the first occurrence of a large 
earthquake within a preexisting GPS network. 

The 1987 earthquake sequence is initially modeled as uniform offsets along 
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rectangular dislocations in an elastic half-space. The best-fit model to the GPS 
observations requires 130 cm right-lateral slip along the northwest-trending Superstition 
Hills fault and 30 cm left-lateral motion along the conjugate northeast-trending Elmore 
Ranch fault. The slip distribution along each fault is investigated by partitioning the 
rupture planes into 50 sub-elements and utilizing Singular Value Decomposition to 
estimate the slip along each sub-fault. Measurements of surface offset are used to 
constrain slip on the shallow elements of the fault plane. The estimated slip distribution 
along the Elmore Ranch fault is fairly uniform. Slip along the Superstition Hills fault 
appears to be concentrated to the northwest and the southeast, with a displacement drop 
near the faults midsection. The estimated geodetic moments are 9.4 x 10 25 dyne-cm and 
2.3 x 10 25 dyne-cm for the Superstition Hills and Elmore Ranch faults, respectively. 
These moments are consistent with those obtained from teleseismic data, suggesting that 
postseismic slip along the Superstition Hills fault was concentrated near the surface. 

In addition, the 1986-1988 GPS observations suggest non-seismic movements across 
the Imperial Valley of up to 5.9 cm/yr. These secular displacements are attributed to 
plate-boundary deformation due to the relative motion between the North American and 
Pacific plates. The observed rate is probably an overestimate, however, as it is heavily 
influenced by unmodeled seismic effects and measurement error. Regardless, the observed 
seismic and secular deformations clearly emphasize the importance of future GPS study 
in the Imperial Valley. 
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Table 1. Imperial Valley GPS Campaign Summary 


Year 

Month 

Days 

Stations 

Organization 

1086 

May /June 

20 

42 

NGS 

1088 

February /March 

0 

10 

UNAVCO 

1088 

March/April 

6 

21 

NGS 

1000 

April 

1 

3 

UNAVCO/RCFC 


NGS - National Geodetic Survey 
UNAVCO - University Navstar Consortium 
RCFC - Riverside County Flood Control 


Table 2. Parameters for Modeled Dislocations 


Superstition Hills 

Elmore Ranch 

Length (km) 

25 

25 

Width (km) 

10 

10 

Strike 

N50* W 

N40*E 

Dip 

00 

00 

Depth (km) 

0 

0 

Latitude ( * N) 

32.0560 

33.1078 

Longitude ( * E) 

-115.7431 

-115.7505 


Latitude and Longitude are coordinates 
at top center of dislocation 
Depth is depth to top of fault 


Table 3. Inverse Models of Seismic Slip 


Model 

Fault 

Sub-faults 
Strike Dip 

Slip 

(cm) 

Moment 
(xlO 25 dyne-cm) 

Model 1 

SH 

1 

1 

100. ± 13. 

7.8 

Model 1 

ER 

1 

1 

-45. ± 10. 

3.4 

Model 2 

SH 

1 

1 

130. ± 8. 

0.4 

Model 2 

ER 

1 

1 

-30. ± 10. 

2.3 

Model 3a 

SH 

10 

5 


0.0 

Model 3a 

ER 

10 

5 


5.0 

Model 3b 

SH 

10 

5 


8.4 

Model 3b 

ER 

10 

5 


7.0 

Model 3c 

SH 

10 

5 


6.2 

Model 3c 

ER 

10 

5 


3.0 

Model 3d 

SH 

10 

5 


0.2 

Model 3d 

ER 

10 

5 


4.0 


SH - Superstition Hills fault 
ER - Elmore Ranch fault 




Table 4. Moment Comparison 


Method Moment (xlO 86 dyne-cm) Reference 



SH 

ER 

Total 

Ratio 


GPS (Model 2) 

9.4 

2.3 

11.7 

4.1 

This Study 

Tdeseismic 

7.2 

1.4 

8.6 

5.1 

Dzicwoncki ct al. [1888] 

Teleseismic 

10. 

2.8 

12. 

4.3 

Sipkin [1989] 

Tdeseismic 

10.8 

2.7 

13.5 

4.0 

Bent ct a 1 . [1989] 

Tdeseismic 

8. 




Hwang ct al. [1990] 

Strong Motion 

5.2 




Wald ct al . [1989] 

Strong Motion 

1.8 




FrankeU and Wcnncrbcrg [1989] 

Pinyon Flat (Planar-A) 

8.7 

0.8 

4.3 

4.6 

Agncw and Wyatt [1989] 

EDM 

9.8 




Lisowekt and Savage [1988] 


SH - Superstition Hills fault 
ER - Elmore Ranch fault 




Figure Captions 


Figure Is Seismicity and major fault systems of the Imperial Valley. The Brawley 
Seismic Zone is the region of anomalously high activity between the northern 
Imperial and southern San Andreas faults. The valley represents a transition zone 
between crustal spreading in the Gulf of California to the south and right-lateral 
transform motion along the San Andreas system. 

Figure 2: Seismicity and surface faulting associated with the November 24, 1987 
Superstition Hills earthquake sequence. A Mg 6.2 event occurred along a northeast- 
trending structure (referred to here as the Elmore Ranch fault) and was followed by 
12 hours with a M s 6.6 event along the northwest-trending Superstition Hills fault. 
The focal mechanisms, aftershock distribution, and surface oflset measurements are 
consistent with left-lateral strike-slip motion along the Elmore Ranch fault and 
right-lateral strike-slip motion along the Superstition Hills fault. A significant 
amount of postseismic slip was observed along the surface trace of the Superstition 
Hills fault, while activity essentially ceased on the Elmore Ranch fault after the M s 
6.6 event. The shaded strips along each fault indicate the geometrical extent of the 
dislocations used to model the geodetic displacements. 

Figure 3: a) Imperial Valley and other GPS stations surveyed in 1986, 1988, and/or 
1990. b) Central-southern Imperial Valley GPS stations surveyed in 1986, 1988, 
and/or 1990. 

Figure 4: Imperial Valley GPS station displacements between 1986 and 1988. All 
movements are relative to station OCTI. The observed displacements are attributed 
to the 1987 earthquake sequence, secular plate-boundary deformation across the 
Imperial Valley, and measurement error. Movements near the 1987 rupture zone 


approach 0.5 m. Error ellipses are determined by multiplying the formal errors by a 
variance factor, determined so the average error scales as 1 ppm (parts per million). 
The east-trending uncertainty is about 4 times larger than the north-trending 
uncertainty. The anomalous southwest-trending apparent movements for those 
stations to the southeast are attributed to measurement error. 

Figure 5: Imperial Valley GPS station displacements between 1686 and 1690. The 
displacement at MOUN relative to BLAC is shown by the dashed arrow. For 
consistency with the 1986-1988 observations the estimated displacement at MOUN 
relative to OCTI is calculated by subtracting the estimated secular velocity of BLAC 
relative to OCTI obtained from conventional geodetic measurements (see Figure 9). 
Station MOUN was reset between surveys and site inspection during 1990 revealed a 
positioning error of about 15 cm. 

Figure 6: Imperial Valley GPS station displacements between February /March 1988 
and March/April 1988. The observed movements indicate the magnitude of errors 
due to the 1988 survey. The vector scale is twice that of Figures 4 and 5. The 
displacement at KANE could represent postseismic deformation from the 1987 
earthquake sequence. 

Figure 7 : The best-fit solution to Model 1 (3-station inversion). The solid arrows 

indicate the observed displacements, while the dashed arrows represent the computed 

displacement based on Model 1. The shaded region indicates where horizontal 

displacements ar^ greater than 4 cm. 

* 

Figure 6: The north and east-trending GPS displacement components for the 1986- 
1988 interval. All distances are computed relative to OCTI on a cross section 
trending N50*E, perpendicular to the plate motion direction. The effects of the 



1987 Superstition Hills earthquake sequence based on Model 1 are removed. Open 
circles indicate stations -where the seismic correction is greater than 4 cm (for that 
component). The 8.1 cm north-trending offset between stations on opposite sides of 
the valley is equivalent to 5.9 cm/yr displacement oriented N40*W. The large 
scatter for the east-trending component is presumably due to large measurement 
errors in the 1986 survey. 

Figure 0: Imperial Valley EDM station velocities computed between 1972 and 1987 
(solid arrows). The movements are largely attributed to secular deformation due to 
the relative motion between the North American and Pacific plates. The 
displacements are modeled (dashed arrows) by 3.45 cm/yr displacement across the 
valley, with right-lateral simple shear oriented N40* W occurring in a transition 
zones 50 km wide north of the Imperial Fault and 20 km wide to the south (shaded). 
The secular deformation is subtracted from the GPS displacements shown in Figures 
4 and 5. 

Figure 10s The residuals (observed minus calculated) for the best-fit solution to 
Model 2. The large residual components in the east-trending direction are suggestive 
of measurement error. The residual at MOUN is likely attributed to the reset 
benchmark between surveys. The unusually large discrepancy at L589 (and KANE) 
suggest additional seismic deformation not accounted for by the simple uniform slip 
parameterization considered for Model 2. 

Figure 11* Observed surface slip along the Superstition Hills and Elmore Ranch 
faults (dashed lines). The Superstition Hills offsets were measured on January 25 and 
26, 1988 about 1 month before the GPS observations. Decaying afterslip is recorded 
up to nearly 1 year after the earthquake sequence. The Elmore Ranch measurements 
are the cumulative slip from multiple surface breaks across a 10 km wide zone, with 


no recorded postseismic offset after the earthquake. The discrete approximation to 
the surface slip used to constrain the uppermost sub-fault elements in Models 3a and 
3b is shown by the solid lines. 

Figure 12: The geodetic moment, standard deviation of sub-fault slip (instability), 
and solution RMS calculated for different singular values (k). a) Shown are solutions 
constrained by surface slip measurements (20<k<30). b) Shown are solutions 
unconstrained by surface slip measurements (1 <k <10). 

Figure 13: The horizontal slip distribution calculated independently for the 
Superstition Hills and Elmore Ranch faults based on Model 2. The shaded region 
indicates where the horizontal deformation is greater than 2 cm for the Superstition 
Hills fault and 0.5 cm for the Elmore Ranch fault. The scale for the Elmore Ranch 
event is altered to account for the smaller dislocation. The deformation pattern is 
almost identical between the two faults, although the displacement magnitudes are 
larger for the Superstition Hills event. This illustrates the difficulty in using the 
GPS measurements to resolve slip between the two faults. 

Figure 14: Slip distribution along the Superstition Hills and Elmore Ranch faults 
computed from Singular Value Decomposition. Each fault is partitioned into the 50 
sub-elements indicated by the grid spacing, a) Solutions for k — 23 and k — 27 
constrained by measurements of surface offset, b) Solutions for k *3 and k*7 
unconstrained by surface slip measurements. 

Figure 15: Known and/or potential conjugate/cross-fault seismic episodes in the 
Imperial Valley since 1900. Seismic release on (left-lateral) northeast-trending 
structures was observed in 1979, 1981, and 1987. Earthquake pairs or 
mainshock/aftershock sequences suggestive of conjugate faulting were observed in 



1015, 1927, 1940. This suggests conjugate/cross-fault interaction is typical for the 
Imperial Valley. 

Figure 18: Strike-shear (right-lateral positive) and normal stress change (dilatation 
positive) induced on the Superstition Hills fault due to a 30 cm left-lateral Elmore 
Ranch dislocation (Model 2). Also shown is the Coulomb failure stress change, 
where positive values indicate an increased potential for rupture (earthquake 
inducing stress). The northwest third of the 1987 Superstition Hills rupture plane 
underwent a stress change tending it toward failure with the maximum change 
calculated in the epicentral region near the intersection with the Elmore Ranch fault. 
There was no rupture to the northwest where the Coulomb failure stress was 
negative (reduced earthquake potential). The magnitude of the stress change (in 
bars) is in the range of typical earthquake stress drops. 

Figure 17: Schematic of potential earthquake failure processes in the Imperial 
Valley, a) Earthquake failure occurs after some critical stress is reached, b) 
Earthquake failure occurs following a time dependent delay after critical stress 'is 
exceeded. 
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